ut 


fo-tae?  2S1 

UNCLASSIFIED 


NUHERICM.  SCHEHES  FOR  THE  ESTIIMTIOH  OF  FIMCTIONM. 
PRRRHETERS  IN  DISTRIBU. .  <U>  BROHN  UNIV  PROVIDENCE  RI 
LEFSCHETZ  CENTER  FOR  DVNBHICRL  SVSTE. . 

H  T  BANKS  ET  RL.  OCT  8S  LCDS-SS-27 


F/B  8/8 


AD-A167  261 


a3 

o 


IKCUNITV  CLAtSiriCATlON  Of  THIS  **01  fVkan  0«« 


REPORT  DOCUMENTATION  PAGE 

'  -*wb*.'TR-  86-  0  260 


'I.  OOVT  ACCISSIOM  MO 


A'  Title  rand  Sua«ii*j 

Numerical  Schemes  for  the  Estimation  of  Functiona 
Parameters  in  Distributed  Models  for  Mixing  MechaiJ 
in  Lake  and  Sea  Sediment  Cores 


AUTMOMfaJ 


H.T.  Banks  and  I.G.  Rosen 


mcmeommino  omsamization  name  ano  aoonsss 

Lefschetz  Center  for  Dynamcial  Systems 
Division  of  Applied  Mathematics 
Brovm  University,  Providence,  RI  02912 


.  controlling  office  name  and  aooncss 

AF0SR//s/rv\  vrtaYc 

Bolling  Air  Force  Base 
Washington,  DC  20332 


t.  monitoring  agency  name  a  AOORCSVK  dlMaml  <Mi  CatimtlliiA  OMIca) 


RBAO  INSTRUCTIONS 
BBrORE  COMPLETING  PORN 


1.  RBCINlCNT't  catalog  NUMGER 


f-  tyre  of  rbrort  a  rcrioo  COVCREO 

isms _ 


$■  CoNtRAcV  OH  grant  NUMGERfai 

oJ^S 

AFOSR-84 


10.  RROGRAM  CLEMENT.  RROJCCT.  TASK 
AREA  A  RORK  UNIT  NUMOCRS 


*1.  RCRORT  DATE 

October  1985 


l>.  NUMGER  OF  RAGES 


It.  security  class.  r*«  thi* 
unclassified 


Is*.  DCCL  AttlFICATION/  OORNGRAOINC 
SCnEOULE 


It  OISTRIGUTION  statement  C*!  Ml*  R*r*/0 


Approved  for  public  release:  distribution  unlimited 


IT-  DiSTRIGuTION  statement  CA*  (Aairaci  atilaraU  In  GlacA  20.  If  Alffaranl  Raai  RaRaRj 


10.  surrlEmEntary  notes 


DTIC 

ELECTE 

tR!T-trt 


IS.  KEY  RORDS  rconiinu*  an  rarara*  alOa  If  nacaaaarr  anU  lUanlllF  Or  Alac*  niaiAar^ 


20.  AGSTR  ACT  rCanlInu*  an  raaaaaa  alR*  If  nacaaaarr  R<F  IRaniffT  *r  Afac*  wiAarf 

- We  consider  distributed  parameter  models  for  vertical  mixing  in  lake  and 

sea  sediment  cores.  Finite  dlrnenslonal  approximation  schemes  are  developed 
for  the  solution  of  associated  Inverse  problems.  The  schemes  permit  one  to 
estimate  temporally  and  spatially  varying  functional  parameters  which  appear 


I® 


DO  1473  COITION  OF  I  NOV  tS  IS  OOBOLBTB 


S/N  0)02'  Lr-014-«A0I 


Unclassified 


I  • 


l; 

i.' 


d 


SBCumrv  CLASStFICATION  OF  THIS  RAGE  (WItm  Oaf*  EnfaraR) 


I 

s 

e| 

S;: 

#  I 


AFOSR-TR.  8G-  0  260 


NUMf.KlCAL  SCHEMAS  FOR  THE  ES 1  IMAM ilN'  OF 
KL'NrTU)NAL  PAR.\ML1KKS  IN  DISTR  1  Rtl  1  EU  MOFKI.S 
FOl'^  MIXINC  MECHANISMS  IN  TAKE  AND  SEA  S!:DlMi:N;  COSES 

hv 

H.  T.  baiU-s  and  C. 


liaerlcal  Scheaes  for  the  Estimation  of  Functional  Parameters  In 


Distributed  Hodela  for  Mixing  Heehanta—  in  Lake  and  Sea  Sedleent  Corea. ^ 


H.T.  Banks* 

Lefschetz  Center  for  Dynamical  Systems 
Division  of  Applied  Mathematics 
Brown  University 
Providence,  RI  02912 


and 

_  _  _  ** 

l.G.  Rosen 

Department  of  Mathematics 
University  of  Southern  California 
Los  Angeles,  CA  90089 


ABSTRACT 


We  consider  distributed  parameter  models  for  vertical  mixing  in  lake  and 
sea  sediment  cores.  Finite  dimensional  approximation  schemes  are  developed 
for  the  solution  of  associated  inverse  problems.  The  schemes  permit  one  to 
estimate  temporally  and  spatially  varying  functional  parameters  which  appear 
in  the  parabolic  partial  differential  equations  and  boundary  conditions 


constituting  the  models.  Theoretical  convergence  results  are  established. 
Numerical  findings  are  presented  which  demonstrate  the  potential  of  the 
methods.  An  example  involving  the  identification  of  a  depth  dependent  mixing 


parameter  based  upon  volcanic  ash  data  from  the  North  Atlantic  is  Included. 
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1.  IWntODOCTKMI 

In  this  paper  we  develop  numerical  approximation  methods  for  the  solution 
of  inverse  problems  arising  in  the  modeling  of  mixing  mechanislms  evidenced  in 
lake  and  sea  sediment  core  samples.  In  particular  we  develop  methods  for  the 
estimation  of  temporally  and  spatially  varying  parameters  in  models  involving 
distributed  (partial  differential  equation)  systems.  Our  work  has  been 
motivated  primarily  by  the  recent  Interest  in  models  involving  depth  dependent 
mixing  rates  and  time  dependent  sedimentation  rates.  These  ideas  will  be 
discussed  in  detail  in  the  next  section. 

The  approach  to  be  described  below  represents  a  significant  Improvement 
over  an  earlier  treatment  of  the  problem  by  the  present  authors  (see  16J). 
Indeed,  based  upon  a  weak  or  variational  formulation  of  the  partial 
differential  equation  and  boundary  conditions  which  compose  the  model  that  is 
the  focus  of  our  study  (as  opposed  to  using  a  semigroun  theoretic  approach), 
the  theory  developed  here  leads  to  schemes  which  permit  the  estimation  of  both 
temporarily  and  spatially  varying  parameters  under  rather  mild  assumptions  on 
the  set  of  admissible  parameters.  No  prior  assumptions  need  be  made  about  the 
"shape"  of  the  functional  parameters  being  estimated. 

Another  feature  of  the  estimation  problems  that  are  of  interest  to  us 
here  which  must  be  considered  when  approximation  schemes  are  being  developed 
is  that  they  generally  Involve  observations  at  the  boundary.  This  noses  a 
significant  mathematical  problem  since  many  mixing  models  involve  parabolic 
systems  which  are  often  formulated  in  the  state  space  where  point 
evaluation  is  undefined.  Under  rather  mild  regularity  assumptions  we  are  able 
to  argue  that  our  state  approximations  converge  in  the  stronger  and,  via 

the  Sobolev  embedding  theorem,  C  topologies.  Consequently,  our  schemes  are 


applicable  to  inverse  problems  in  which  the  fit  to  data  criterion  Involves 


point  or  discrete  as  well  as  distributed  observations  In  the  spatial  variable. 

The  general  approach  we  take  here  has  been  used  extenslvelv  to  develop 
approximation  methods  for  the  solution  of  Inverse  problems  Involving 
distributed  systems  arising  In  a  wide  varletv  of  application  areas  (e.g. 
population  dispersal,  physiology,  large  flexible  spacecraft,  seismic  analvsls, 
diffusion  through  porous  media,  etc.).  A  brief  description  of  some  of  the 
problems  and  the  associated  approximation  schemes  together  with  a  survey  of 
the  literature  and  a  more  comprehensive  blbllographv  can  be  found  In  {IJ. 

We  provide  a  brief  outline  of  the  remainder  of  the  paper.  In  Section  2 
the  sediment  core  mixing  problem  Is  described  and  previous  modeling  efforts 
are  discussed.  A  particular  model  Involving  a  dlf fuslon/advectlon  equation  Is 
developed  In  detail  and  associated  Inverse  problems  are  posed.  The 
Identification  problem  Is  formally  stated  and  the  weak  formulation  of  the 
partial  differential  equation  together  with  existence,  uniqueness  and 
regularity  results  for  solutions  are  discussed  In  Section  3.  In  Section  4  we 
describe  the  approximation  schemes  and  establish  convergence  results. 

Examples  and  our  numerical  findings  are  given  in  Section  5.  Our  primary 
objective  in  this  paper  Is  to  present  our  theoretical  results.  The  numerical 
results  discussed  In  Section  5  are  preliminary.  A  more  complete  numerical 
study  will  be  described  elsewhere. 

The  notation  we  employ  Is  standard  throughout.  We  denote  by  H'^(a,b)  the 

usual  Sobolev  spaces  of  real  valued  functions  ^  defined  on  the  Interval  (a,b) 

k~l  v» 

whose  (k-l)st  derivative,  D  Is  absolutely  continuous  and  whose  k^" 

k  0 

derivative,  D  ^  is  In  L2(a,b)  ■  H  (a,b).  The  standard  Sobolev  Inner  products 
and  norms  are  denoted  by  respectively.  For  *  H'^(a,b)  a 

function  V  ;  (t^,!^)  -►  is  said  to  be  an  element  In  H®(tQ,tj;  h'^)  If 
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4  V 

We  shall  say  v  c  H''(tQ,t^;H  )  if  v  posseses  1  stront?  derivatives  with 
v^^^  t  L2(tQ,t^;H'"),  0  <  1  <  j.  Defining 


j  t 


i  k  -  ^  ^  /t 
'•  1-0  0 


1  .m2  ...1/2 


)l,dt}' 


*  If  ^ 

for  V  i  HJ(tQ,tj;H  )  we  have  that  the  spaces  {H  (t^.t^H  ),  ar® 

complete.  For  j-k,  we  abbreviate  notation  for  the  norm  by  writing 

1*1  ^1*1 
'  'k  '  'k,k* 

Finally  for  V  a  Hilbert  space  which  Is  densely  and  compactly  embedded  In 
and  V  c  H^(tQ,tj;V)  we  denote  the  weak  derivative  of  v,  an  element  In 
H®( tQ, t £ ; V ' ) ,  by  V  where  V  Is  the  dual  space  of  V.  For  i  t  V  we  have 


<v(c))(<»  -  <^,v(t)X 


where  the  Inner  product  In  the  above  expression  Is  understood  to  be  the 
extension  of  the  Inner  product  to  all  of  V. 


2.  Inverse  ProblcM  In  Lake  and  Sea  Sediment  Analysts 

Sediment  formation  in  lakes  and  deep  seas  is  of  threat  importance  to 
geophysical  scientists  who  use  core  samples  of  sediment  in  their 
investigations  of  the  history  (e.g.  palaeoc lima tic  changes)  of  the  earth. 
Unfortunately,  the  stratigraphic  records  contained  in  these  core  samples  have 
been  subjected  to  perturbations  since  ocean  and  lake  floors  are  in  general  not 
quiescent.  Two  general  types  of  redistribution  of  sediment  are  often 
significant:  (1)  gross  lateral  transport  via  ocean  and  lake  bottom  currents, 
for  example,  though  a  continuous  winnowing  of  bottom  currents  or  through 
episodic  currents  such  as  turbidity  currents;  (ii)  the  mixing  activities  of 
benthic  organisms  near  (on  the  order  of  2-40  cm.)  the  sediment-water 
Interface.  This  biological  mixing  of  sediments  by  organisms  (which  leads  to 
an  interesting  class  of  inverse  problems  in  the  analysis  of  sediments)  is 
called  bioturbation  124J,  (7J  and  takes  place  in  sediment  layers  in  bodies  of 
water  (lakes,  estuaries,  the  deep  oceans)  in  which  bottom  water  is  not 
substantially  depleted  of  oxygen.  Bioturbation  is  effected  bv  different  kinds 
of  organisms  such  as  clams,  worms,  Crustacea,  echlnoderms,  etc.,  and  the 
mixing  activities  consist  primarily  of  burrowing  (e.g.,  for  safety)  and 
ingestion  -  excretion  reworking  of  the  sediment  for  its  edible  organic 
matter.  Through  the  use  of  tracers  from  dated  events  (e.g.,  nlutonlum-from 
atmospheric  fallout  from  nuclear  exposlons,  and  mlcrotektites  -  tiny  drops  of 
sculptured  glass  resulting  from  cosmic  events),  it  can  be  determined  that  the 
biological  mixing  of  abyssal  sediments  is  quantitatively  significant  and  takes 
place  on  a  relatively  short  (in  regard  to  geologic  records)  time  scale  (10-20 
years).  Furthermore,  there  seems  to  be  little  correlation  between 
bioturbation  mixing  rates  (which  are  highly  variable)  and  sediment  type  or 
sediment  accumulation  rates.  Howeyer,  the  degree  of  bioturbation  and  the 


depth  of  the  region  in  which  it  occurs  are  related  to  the  types  of  organisms 
inhabiting  a  particular  area. 

Since  bioturbation  plays  such  a  fundamental  role  in  the  alteration  of 
geologic  records,  it  is  not  surprising  that  geochemists,  geologists,  and 
geophysicists  have  in  recent  years  attempted  to  understand  the  effects  of 
bioturbation  well  enough  so  as  to  enable  one  to  properly  interpret  the 
Information  contained  in  core  samples,  thereby  sharpening  the  details  in  these 
geologic  records.  A  number  of  Increasingly  sophisticated  mathematical  models 
along  with  related  "inverse  problems"  can  be  found  in  the  literature  111)  112| 
[13J  115J  [161  [17]  [18]  [19].  These  models  typically  involve  some  tvpe  of 
region  or  chamber  (ranging  from  a  simple  well-mixed  chamber  to  one  in  which 
mixing  rates  are  depth  dependent)  in  which  mixing  and  advection  or  convective 
flow  interact  to  vertically  redistribute  sediment  particulate  matter,  volcanic 
ash,  microtektltes,  radioactive  tracers  or  other  substances  from  episodic  and 
nonepisodic  events. 

One  model  which,  along  with  its  variations,  has  en.loyed  rather  widespread 
usage  involves  the  assumption  that  one  has  a  vertically  moving  chamber 
(assumed  uniform  in  horizontal  directions)  in  which  mixing  and  advectlve  flow 
of  material  takes  place  and  is  described  by  one-dimensional  (depth)  transport 
equations.  Depth  in  the  chamber  is  represented  by  coordinates  x,  0  <  x  <  I 
and  the  chamber  (and  hence  coordinate  system)  is  assumed  to  be  moving  upward 
with  a  velocity  ■  l^(t)  (corresponding  to  sedimentation  rate  or,  equivalently 
in  this  case,  sediment  layer  buildup)  so  that  it  is  always  located  in  the 
top  I  cm.  of  the  sediment  as  depicted  in  the  figure  below.  Thus  x  *  0  is 
always  at  the  water-sediment 


water 


1  advective 

material 

flux 

sediment 

V  =  (/(t) 

i 


0 


I 


Figure  2.1 


interface  and  the  bottom  of  the  chamber  at  x  =  i  is  located  at  that  denth 
beyond  which  (it  is  assumed)  no  further  changes  (i.e.,  no  bioturbation)  in  the 
historical  records  occur.  The  resulting  configuration  with  upward  velocity 
of  the  chamber  can  be  equivalently  modeled  by  the  assumption  of  a  fixed 
coordinate  system  for  the  chamber  with  an  advective/conyective  flow  of 
material  downward  through  the  chamber  with  velocity  .  Use  of  the  model  by 
numerous  Investigators  (e.g.  see  (18J,  119J)  strongly  suggests  that  permitting 
a  time-varying  sedimentation  rate  in  such  models  is  Important. 

If  u  =  u(t,x)  is  the  concentration  of  material  (e.g.,  shards  of  ash, 
radioactive  tracer,  etc.)  with  whose  movement  one  is  concerned  and  i  ■  i(t,x) 
is  the  material  flux  at  time  t  and  position  x  in  the  chamber,  material 
conservation  is  represented  by  the  classical  mass  balance  or  continuity 
equation 

<2-‘)  If  ♦ 

where  X  is  a  decay  constant  for  the  material  (X  *  0  if  one  is  dealing  with  a 
conservative  tracer  such  as  mlcrotektites).  Of  course,  the  Important  aspect 
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of  any  such  model  is  the  assumption  one  makes  regarding  the  material  flux  i, 
which  here  we  assume  consists  of  a  mlxlnt?  component  and  an  advecttve 
component.  In  this  case  one  is  justified  in  assuming  that  the  bloturbation 
(burrowing,  ingestion,  etc.)  takes  place  over  a  very  short  time  scale 
(essentially  Instantaneous)  in  the  chamber  and  hence  perhaps  can  be 
represented  by  a  dlf fuslve-llke  flux  component.  The  advective  flux  is  given 
by  t'(t)u  and  if  one  assumes  a  Fickian  flux  for  the  bloturbation  with  depth 
dependent  "bioturbatlon’’  coefficient  P=  P(x),  one  obtains 

(2.2)  j  *  -  ^^§7  +  l/u  . 

The  assumption  that  0  is  a  function  of  depth  is  motivated  bv  ones 
expectation  that  the  rate  of  mixing  is  generally  higher  near  the  well- 
oxygenated,  densely  populated  surface  of  the  sediment  mixing  layer;  this 
expectation  appears  to  be  corroborated  by  experimental  findings  lllj,  |12J, 
113J,  [19J.  A  more  fundamental  question  as  to  whether  the  biological 
reworking  of  sediment  is  mechanistically  analogous  to  molecular  diffusion  (and 
hence  is  consistent  with  the  Fickian  flux  assumption)  is  not  so  readilv 
answered.  A  strict  analogy  would  neceasltate  the  existence  of  an  abundance  of 
organisms,  randomly  placed  in  the  chamber,  mixing  the  materials  in  a  manner  so 
as  to  produce  a  material  flux  proportional  to  concentration  gradients.  While 
this  is  not  a  very  likely  description  of  the  mechanisms  of  biogenic  mixing, 
one  might  still  have  a  plausible  quantitative  analogy  with  diffusion  if  the 
mixing  rate  is  rapid  and  sediment  samples  which  Involye  a  large  number  of 
Independent  transport  events  of  variable  duration  are  chosen. 

For  boundary  conditions  at  the  upper  boundary  (x  =*  0)  of  the  chamber  one 
has  the  flux  condition  j(t,0)  *■  G(t)  when  G  is  a  possibly  unknown  input,  while 
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the  total  flux  at  the  lower  boundary  x  =  £  is  via  advectlve  loss  throu<?h  the 
bottom  of  the  mixing  zone  and  thus  j(t,£)  =  l/(t)u(t,£).  Using  the 
constitutive  relationship  for  j  from  (2.2)  in  these  boundary  conditions  and  in 
the  equation  (2.1),  one  obtains  the  model 

(2.3)  I7  =  (P(x)  I7)  -  l'(t)  1“  -  Xu,  0  <  X  <  £,  t  >  0, 

(2.4)  -  P(0)  1^  (t.O)  +  l/(t)u(t,0)  =  G(t) 

(2.5)  -  P(£)  (t,£)  =  0 

(2.6)  u(0,x)  =  $(x) , 

where  4  is  the  initial  distribution  of  material  in  the  chamber. 

The  appropriate  initial  data  assumption  is  closely  related  to  the 

assumption  one  makes  about  the  input  flux  G.  To  Illustrate  possibilities,  we 

can  consider  several  specific  situations  that  arise  in  geological 

investigations.  A  strong  argument  for  steady-state  input  flux  (G(t)  = 

constant)  can  be  made  in  the  case  one  is  investigating  a  tracer  such  as  lead 
210 

(  Pb)  which  exhibits  a  rather  steady  production  rate  from  atmospheric  (decay 
of  gaseous  radon  -  222)  and  oceanic  (decay  of  radon-226)  sources.  For  tracers 
such  as  plutonium  (^  ’  ^Pu)  and  cesium  (*^  'Cs)  which  result  from  atmospheric 

nuclear  weapons  testing,  time  dependent  flux  is  more  anpropriate,  and  in  some 
cases,  for  appropriately  chosen  initial  times,  the  assumption  that  ♦  vanishes 
is  appropriate.  In  general  though,  in  both  cases  one  must  estimate  the 
initial  distribution  either  as  a  part  of  the  overall  inverse  problem  or, 
through  an  a-priori  procedure  using  directly  earlier  (i.e.  deeper) 


concentration  profiles  In  the  sediment  core.  Finally,  in  the  case  of  truly 
episodic  events  (ash  shards  from  volcanic  eruptions,  microtektltes  of  cosmic 
origin)  an  Impulse  input  is  most  appropriate.  This  can  be  effectively  modeled 
by  choosing  an  impluse  like  initial  function  ♦  in  (2.6)  and  taking 
G  »  0  in  (2.4).  The  magnitude  of  this  Impulse  can  sometimes  be  rather  easily 
estimated  directly  from  knowledge  of  the  total  material  content  of  the  sample. 

In  any  case,  to  understand  the  effects  of  bloturbatlon  on  the 
distribution  of  material  concentrations  in  core  samples,  it  is  sufficient  to 
have  values  for  the  parameters  V,  V,  X,  and  A,  and,  of  course,  to  know  that 
use  of  these  parameter  values  in  the  model  gives  one  an  accurate  quantitative 
description  of  concentrations  found  in  core  samples.  It  can  be  expected  that 


these  parameter  values  will  vary  depending  on  the  core  sample  and  the  material 
under  investigation.  Hence  one  would  like  to  have  a  procedure  whereby  given 
data  from  a  specific  core  sample,  one  can,  with  some  confidence,  determine  the 
"correct"  parameter  values.  In  regard  to  this  Inverse  procedure,  for  the 
model  above  we  note  that  concentrations  in  the  historical  layers  (where  time 
in  klloyears  can  usually  be  related  to  centimeters  of  thickness  of  core  sample 
e.g. ,  see  111])  represent  concentrations  at  various  times  at  the  bottom  of  the 
mixing  chamber  (i.e.  at  x  =  A).  Hence  data  for  the  process  may  be  given 
by  Z(C)  where  Z(C)  denotes  the  observed  concentration  of  tracer  material  at  a 
height  of  C  cm.  above  the  position  in  the  core  designated  as  time  t  *  0.  In 
this  event  a  typical  inverse  problem  might  be  stated:  Given 

observations  Z(C^),  1  =  1,2,.. .,<  at  core  locations  1  «  1,2 . k  find, 

among  some  class  Q  of  admissible  parameters,  parameters  q  >  (  P,  l/,X,A)  that 


minimize 


(2.7) 


J(q)  -  I  |Z(C.)  -  u(T(c^;l'),<i)r 
i=l 


^  •  •  •  •  a  •  fc**  .*• 


O  m.**  •*.  . 
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where  u  is  the  solution  to  (2.3)-(2.6)  corresponding  to  q  and 
T(5;  I')  -  r“^(C)  with  r^^  (t)  -  /J  l/(s)ds. 

In  addition  to,  or  in  place  of  some  of  the  parameters  in  q,  sometimes  it 
is  also  desirable  to  estimate  G  and/or  9  in  the  formulation  above.  Of 
course,  in  some  instances  the  data  from  core  samples  will  not  support  such  an 
inclusive  inverse  procedure. 

The  model  formulated  above  is  based  on  the  assumption  that  the  entire 
chamber  is  available  for  throughput  of  material  and  that  the  sedimentation 
rate  is  the  same  as  the  material  velocity  through  the  chamber  (l.e.  no 
compactif Icatlon  of  sediment  takes  place).  In  many  instances  porosity  effects 
and/or  compactif icatlon  are  important  and  should  be  included  in  the  model.  It 
is  also  sometimes  important  to  distinguish  between  tracer  materials  and 
sediment  particles.  These  concepts  require  modifications  of  the  modeling 
ideas  presented  above. 

We  again  postulate  the  moving  mixing  chamber  but  (for  reasons  that  will 
become  clear  in  the  sequel)  now  we  let  z  denote  the  chamber  coordinates  (z  ■  0 
is  the  water-sediment  interface,  z  ■  £  is  the  bottom  of  the  mixing 
chamber).  The  porosity  ♦  is  the  fraction  of  the  chamber  volume  that  is 
available  for  flow  (throughput)  so  that  1  -  ^  is  the  fraction  that  is 
solid.  We  assume  that  the  porosity  A  >  ^(z)  is  depth  dependent  and  let  p 

s 

be  the  constant  sediment  particle  density  (in  mass  per  unit  length  of 

particulate  matter).  If  we  furthermore  let  l/(t,z)  be  the  sediment  oartlcle 

velocity  with  respect  to  the  z  coordinate  system,  we  may  write  separate  mass 

balance  equations  for  the  sediment  particulate  matter  and  tracer.  Considering 

first  sediment  particles,  we  have  that  the  particle  mass  density  in  the 

chamber  is  given  by  p  (1-^)  and  the  particle  mass  flux  is  given 

s 

by  jg  *  jg(t,z)  ■  Pg(l-^)l/  •  Assuming  conservative  particulate  matter,  we 


V/i 
■-  .V.' 
■•‘.V,  ( 
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obtaln  the  sediment  mass  balance  equation 

(2.8)  (Pg(l-((.))  +  (Pg(l-i)l/)  -  0. 

Note  that  this  conservation  law  Implies  jg  Is  actually  Independent  of  z. 

For  the  tracer,  we  have  similar  considerations  regarding  porosity  excent 
the  flux  Is  more  Involved  since  we  have  bloturbatlon  and  convective  movement 
along  with  decay.  (The  sediment  particles  may  also  be  mixed  (Ingested,  etc) 
but  we  cannot  measure  (observe)  this  ~  our  observations  being  of  tracer 
material.  Hence  the  flux  for  sediment  particles  only  contains  a  convective 
movement  term.)  Let  c  denote  the  mass  of  tracer  per  unit  mass  of  sediment 
particulate  matter.  Then  the  tracer  mass  density  (mass  per  unit  length  of  the 
chamber)  Is  given  by  cp^(l-^)  and  the  tracer  mass  balance  equation  can  be 
written 

(2.9)  1^  (cPgd-^))  +  (j^)  +  X(cp^(l-^))  -  0 

where  Is  the  tracer  mass  flux.  Denoting  the  tracer  velocity  by  1/^  that 
the  tracer  mass  flux  Is  glyen  by  ,1^  •  cp^(l-4>)l/^,  we  may  divide  the  tracer 
mass  flux  Into  components  representing  mixing  (bloturbatlon)  and  convective 
flux 


-  cPg(l-^)(Mj.  -  I/)  +  cpg(l-4.)l/  . 

The  term  ■  cp  (  -  I/)  may  be  regarded  as  a  "pure"  mixing  tracer  flux 

M  8  T 

(l.e.  neglecting  porosity)  and  If  we  make  the  Flcklan  assumption  for  this  flux 


term 


the  resulting  mixing  flux  is  given  by  P"!" 


Hence  we  have 


-  -p,(i-0  P|f  * 

and  the  tracer  mass  balance  equation  (2.9)  can  be  written 

(2.10)  |-  (cp  (!-♦))  +  f-  (P  (l-«)c{/  )  «  |-  (p  (!-♦)  P~)  -  Xp  (l-(>)c. 

dt  S  dz  s  dz  S  dz  S 

The  sediment  mass  balance  equation  (2.8)  can  be  used  to  modify  the  tracer 
equation.  Observing  that 

(CP  (!-♦))  -  c  (p  (!-♦))  +  p  (1-^) 

dt  S  dt  S  S  dt 

.-c|j  (.,(!-♦) V 

-  ‘'If-  I? 

and  using  this  in  (2.10),  we  obtain 

(2.11)  p  (1-^)  |£  +  p  (1-^)  (p  (1-,^)  d|£)  -  Xp  (l-<>)c. 

S  dt  S  dz  dz  S  dz  S 


In  the  case  of  constant  porosity,  this  equation  reduces  to 
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whlch  Is  Che  same  as  (2.3)  if  c  is  Interpreted  as  volumetric  concentration. 

Equation  (2.11)  can  be  put  in  a  form  similar  to  (2.3)  even  when  the 
porosity  Is  not  constant.  We  define  a  new  Independent  variable  x  »  x(z)  bv 


X  »  /q  p^(l-^)dC. 

Assuming  that  p^(l-^)  >  0,  this  defines  an  invertible  mapping  between  the 

chamber  coordinate  variable  z  and  the  new  variable  x  which  is  sometimes  called 

9  3 

the  "total  sediment  particle  accumulation".  We  then  have  -r—  *  p  (l-4i)  -r—  and 

3z  8  3x 

use  of  this  in  (2.11)  yields 


(2.13) 


3c  .  3c  3  /  r  . 

IT  *  “  57  •  57  '  f 


2 

where  E  =  (p  (l-^i))  V  and  (o  =  p  (l-(j))l/  Is  the  particle  mass  flux  (called 
s  s 

jg  above).  We  see  that  equation  (2.13)  has  the  same  form  as  (2.3)  and  note 
Chat  in  this  equation,  c  can  be  expressed  either  as  tracer  mass  per  unit  mass 
of  sediment  particulate  matter  or  in  the  more  traditional  interpretation  as 
tracer  mass  per  unit  length  of  the  sediment  column.  Also  note  that  in  (2.13) 
X  Is  no  longer  slmplv  the  vertical  distance  in  the  column  nor  is  E  simply  the 
mixing  or  bloturbation  coefficient.  Moreover,  in  such  models  where  porosity 
and  compact iflcat ion  are  considered,  the  concent  of  "sedimentation  rate"  is 
more  delicate.  To  elaborate  on  this,  we  may  define  M(t),  the  total  sediment 
mass  at  time  t  by 


M(t)  -  /^^*'^p^(l-«)d5 
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where  h(t)  is  the  sediment  column  depth  measured  from  some  fixed  reference 

coordinate  system.  Note  that  ^  represents  the  rate  at  which  the  water- 

dt 

sediment  interface  Is  changing  relative  to  a  fixed  coordinate  system  and  hence 
is  the  rate  at  which  the  sedimentation  layer  is  increasing  .  Thus,  this  is 
the  "true  sedimentation  rate"  V  (t)  which  in  general  differs  from  the  rate 
(t,0)  at  which  particles  are  "falling"  and  passing  into  the  chamber.  The 
porosity  is  often  assumed  to  have  the  form 

(Kz)  -  (i>Q  -  ♦pe  +  41^ 

dM 

where  >  4,  if  compact if icat ion  takes  place.  Moreover,  since  u)  ■  is 

U  I  at 

the  sediment  particle  mass  flux  already  defined  as  p  ( 1-4)1^  ,  we  see  that 

s 

a)(t)  -  Pg(l-*(h(t)))  1^-  p^(l-*(z))  l'(t,z). 

Then  for  4q  >  we  haye 

“  p  (l-4.(h(t)))  p  (l-^(O))  " 
s  s 

or,  for  a  thick  sediment  column  (so  that  ♦  "  we  have 

(I)  X  u> 

p  (!-♦,)  P  (!-♦.) 


1/  (t) 


l'(t,0) 
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3,  TUB  IDKHTlglCATION  PROBLEM 

In  this  paper  we  will  concentrate  our  efforts  on  the  develooment  of 
numerical  approximation  schemes  for  the  solution  of  the  parameter  estimation 
problem  involving  the  simplified  bloturbation  model  given  by  (2.3)-(2.6). 
Methods  applicable  to  inverse  problems  for  more  elaborate  models  (e.g.  those 
involving  porosity  and/or  compactif ication  effects)  are  currently  under 
investigation  and  our  findings  will  be  discussed  elsewhere. 

In  this  section  we  give  a  precise  formulation  of  the  identification 
problem  and  briefly  outline  existence,  uniqueness  and  regularity  results  for 
solutions  to  initial- boundary  value  problems  of  the  form  (2.3)-(2.6)  which 
will  be  required  in  the  subsequent  development. 

If  we  make  the  change  of  variable  y  *•  x/£  and  allow  for  the  inclusion  of 
a  slnk/source  term  F  =  F(t,x)  in  the  partial  differential  equation  (2,3)  we 
may,  without  loss  of  generality,  consider  the  system  for  states  v  »  v(t,y)  and 
parameters  q  *  (q^  (y),q2^  t),  f  ( t ,y),g( t),  4>(y) , X, i)  given  by 

(3.1)  1^  v(t,y)  -  ^  (q^(y)  |;^  v(t,y))  -  q^Ct)  v(t,y) 

-  Xv(t,y)  +  f(t,y),  t  >  0,  y  fc  (0,1) 

(3.2)  -qj(0)  |^v(t,0)  +  q2(t)v(t,0)  -  g(t),  t  >  0 

(3.3)  -qj(l)  Ip-  v(t,l)  =0,  t  >  0 

(3.4)  v(0,y)  -  <t>(y)t  y  ^  (0,1), 


where  the  original  states  and  parameters  appearing  in  (2.3}-(2.6)  can  he 
recovered  from  the  relations 

t?(x)  «  A^q^(x/1), 

l/(t)  - 

F(t,x)  -  f(t,x/t), 

G(t)  =  lg(t), 

*(x)  -  «(x/l) 

and 

u(t,x)  ■  v(t,x/£). 

The  least  squares  performance  Index  (2.7)  takes  the  form 

2 

(3.5)  J(q;v)  «  E  |Z(C,)  -  v(t(C, ; £q,) , 1) T 
1-1  ^ 

where  v  is  the  solution  to  (3.1)-(3.4)  corresponding  to  q. 

In  what  Is  to  follow  we  make  the  standing  assumption  that  the  slnk/source 
term  f,  the  boundary  flux  g  and  initial  conditions  ♦  are  known  or  have  been 
estimated  a-prlorl.  The  parameter  vector  q  is  assumed  therefore  to  be  of  the 
form  q  »  (q^ (y) ,q2(t) Treating  this  somewhat  scaled  down  version  of 
the  original  identification  problem  will  capture  all  of  the  essential  features 


of  our  approach  and  at  the  same  time  keep  the  presentation  and  discussion  as 
simple  as  possible.  In  addition,  this  formulation  Is  more  than  adeauate  to 
treat  the  identification  problem  corresponding  to  a  set  of  observations 
Involving  volcanic  ash  data  which  has  been  the  primary  focus  of  our  research 
to  date  (see  [6J  and  Example  5.3  below).  The  necessary  modifications  to  our 
theory  so  as  to  allow  the  estimation  of  the  input  terms  f  and  g  and  the 
initial  data  ^  are  relatively  straightforward  and  have  been  discussed  in 
detail  elsewhere  (see  [3],  (4J). 

We  assume  e  H^(0,1)  and  that  there  exists  a  T  >  0,  sufficiently  large 
for  which  the  mappings  t  g(t)  and  t  f(t,*)  are  in  H^(0,T)  and 
H®(0,T;H°(0,1))  respectively. 

Define 

Q.»  CIO, 11  X  H*(0,T)  X  r'  X  r' 

with  the  usual  product  space  topology  and  let  Q^Oj  x02’‘A*l>c()  satisfy 
the  following  hypotheses: 

(HI)  is  a  compact  subset  of  C10,1J  and  there  exist  constants  u  and  v 

such  that  0  <  u  <  q^(y)  <  v,  y  t  10,1J  for  all 

(H2)  Q2  is  a  compact  (with  respect  to  the  topology)  subset  of  cMo,TI 
with  0  <  w  <  q2(t)  <  v  and  |q2(t)|  <  v  for  all  t  €  10, TJ  and  all 
<i2  ^  ^2’ 

(H3)  A  is  a  compact  subset  of  R^  with  0  <  X  <  v  for  all  X  €  A, 
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(H4)  L  is  a  compact  subset  of  with  0  <  p  <  £  <  v  for  all  £  €  L. 

It  is  clear  that  (HI)  -  (H4)  above  imply  that  Q  is  a  compact  subset  of  0. 

Letting  H  »  H®(0,1)  and  V  *  H^<0,1)  endowed  with  the  standard  Sobolev 
inner  products  we  have  the  usual  dense  embedding  V  ^  He  V ' .  For 
each  t  t  10, TJ  and  q  =  (qj,q2,A,L)  ^ ^  define  the  bilinear 
form  L(t;q):VxV+R  by 

(-(t;q)(x.4')  =  <q^DX,D4i>Q  -  q2(t)<  X.Dt(>>Q 

+  ^  <x.’l'>g  +  q2(t)x(i)'l»(i),  X.'l'  €  V 

and  consider  the  weak  form  of  (3.1)-(3.4)  given  by 

(3.6)  <v<t),<|/>Q  +  L(t;q)(v(t),il>)  =  <f(t),i>.>p  +  g(t)i|;(0),  t  >  0,i(-  <='  V 

(3.7)  v(0)  =  4. 

where  v(t)  =  v(t,*)  and  f(t)  ”  f(t,»). 

Hypotheses  (HI)  -  (H4)  imply  that  for  each  t  i  10, T],  q  t  0  and  x*'l'  t  V 

I  L(t;q)(x,'l')  I  <  4v  1x1 

and  that  L(c;q)(* ,•)  is  strongly  elliptic.  It  follows  therefore  (see  [10], 
Theorem  10.3.4,  (211  Theorem  III.  5.A.)  that  it  is  V-coerclve;  that  is,  there 
exists  an  a  y  R  such  that 
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i.  (t;q)(i|<,»l>)  +  oI'I'Iq  >  61'1'lj, 

6  >  0,  for  all  i|>  t  V  and  t  i  [0,TJ.  If  we  rewrite  the  terns  on  the  right 
hand  side  of  (3.6)  as 

<f(t),ii»>Q  +  g(t)\|j(0)  »  <f(t),i|>>^ 

A 

where  f(t)  =  f(t,«)  +  g(t)5(»)  and  6(«)  denotes  the  dlrac  Impulse  at  zero, 

“  0 

then  f  c  H  (0,T;V')  and  standard  results  (see  [14J,  Theorem  III. 1.2)  yield 
that  the  system  (3.6),  (3.7)  admits  a  unique  solution  v  ^  C([0,T1;H)  with 
V  t  H®(0,T;V)  and  v  c  H°(0,T;V'). 

Under  additional  regularity  hypotheses  on  q| ,  q2,  f,  g  and  4'»  smoother 
solutions  to  (3.6),  (3.7)  can  be  obtained  (see  below).  When  these  solutions 
are  sufficiently  smooth  and  qp  q2.  f,  R  and  ^  are  sufficiently  regular,  they 
coincide  with  either  strong  (defined  in  terms  of  evolution  operators,  see 
[231)  or  classical  (see  llOJ)  solutions  to  the  original  initial  boundary  value 
problem  (3.1)-(3.4).  In  light  of  this,  we  formally  define  the  identification 
problem  as 

(  P)  Find  q  ^  Q  which  minimizes  J(q;v)  given  by  (3.5)  where  v  is  the 
solution  to  (3.6),  (3.7)  corresponding  to  q. 

As  is  typically  the  case  with  Rltz-Galerkin  based  methods,  in  order  to 
demonstrate  that  our  approximation  schemes  converge,  additional  regularity  of 
solutions  to  (3.6),  (3.7)  will  be  required.  We  state  the  following  theorem. 


Theorea  3»1.  Suppose  q2  t  C  [0,Tl.  Let  n  be  a  function  in 
C^dO.Tj  X  [0,1  J)  which  satisfies 

(1)  n(0,y)  =  (Ky),  y  «  (0,1) 

(11)  -qj(0)  1^  ii(t,0)  +  q2(t)n(t,0)  -  q(t),  t  >  0 

(ill)  -qj(l)  1^  n(t,l)  -  0,  t  >  0 

and  define 

h(t,y)  -  f(t,y)  +  1^  {qj(y)  n(t,y)}  -  q2(t)  n(t,y)  -Xn(t,v)  -  n(t 

If  h  c  C^(10,TJ  X  [0,1J)  and 

g.i 

^  h(0,y)  -  0,  y  €  (0,1),  ,1  -  0,1, 

3t^ 

then  V,  the  solution  to  (3.6),  (3.7)  is  an  element  in  H^(0,T;  H^(0,1))  and  it 
can  be  Identified  with  functions  in  H^((0,T)  x  (0,1)). 

The  proof  of  Theorem  3.1,  which  has  been  omitted,  can  be  arqued  usinq 
Theorem  10.6.17  in  [lOJ. 

The  conditions  specified  in  the  statement  of  Theorem  3.1,  which  of  course 
are  only  sufficient  conditions,  are  rather  restrictive.  However,  they  can  be 
used  to  derive  specific  and  not  unreasonable  conditions  on  the  initial  and 
boundary  data  (♦,q)  and  parameters  qj,q2  to  insure  reqularlty  of  solutions  to 
(3.6),  (3.7).  Althouqh  the  theorem  above  is  not  the  most  qeneral  one 
possible,  it  suffices,  for  our  purposes  here,  to  simply  demonstrate  the 


f  'j^  irvv^n.'v^-^  L*^'  v^  v.~»  '."» »>■*  ^ 
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existence  of  H^(0,T;H^(0,1))  solutions.  It  is  in  fact  the  case,  that  while 
regularity  of  solutions  is  necessary  to  demonstrate  convergence,  it  has  been 
our  experience  that  it  has  little  or  no  effect  on  the  actual  performance  of 


the  schemes  when  they  are  applied  in  oractlce 


4.  APPROXIMATION 


In  this  section  we  develop  approximation  schemes  to  solve  the  parameter 

estimation  problem  (P).  Fundamental  to  our  approach  is  the  construction  of  a 

sequence  of  finite  dimensional  identification  problems,  each  of  which  has  a 

solution  that  can  be  found  uslnq  conventional  techniques  and  readily  available 

software.  We  demonstrate  that  under  appropriate  hypotheses,  solutions  to  the 

finite  dimensional  problems,  in  some  sense,  approximate  solutions  to  the 

original  infinite  dimensional  identification  problem  (P). 

A  two  stage  approximation  process  is  employed.  We  first  use  a  Rltz- 

Galerkin  approach  to  approximate  the  infinite  dimensional  state  equation  (3.6) 

by  a  sequence  of  finite  dimensional  ordinary  differential  equations.  The  set 

of  admissible  parameters  Q,  a  subset  of  the  infinite  dimensional  function 

space  0.  is  then  discretized.  The  result  is  a  sequence  of  optimization 

problems  involving  the  minimization  of  a  least  squares  performance  index  over 

a  compact  subset  of  Euclidean  space  subject  to  finite  dimensional  constraints. 

For  each  N  =  1,2,...  let  be  a  finite  dimensional  subspace  of  H 

N 

satisfying  V”  V.  Suppose  V”  =  span  a’ld  define  P  :H  V  to  be  the 

orthogonal  projection  of  H  onto  with  respect  to  the  H  (i.e.,  <»,*>^)  inner 

product.  We  make  the  following  hypothesis  about  the  approximating  properties 
N 

of  the  subspaces  V  . 

(H5)  For  each  t);  t  h\o,1),  +  0  as  N  •  with  <  kl<l)|j  for 

some  constant  k  independent  of  N  and  <ti. 

The  usual  spaces  of  linear  or  cubic  B-spllne  functions,  among  others,  are 
known  to  satisfy  Hypothesis  (H5)  (see  (20),  (22 J). 
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The  Galerkin  equations  in  corresponding  to  the  system  (3.6),  (3.7)  are 
given  by 

(4.1)  <  v^(t),i(/^  >Q  +  (.(t;q)(v^( t),i;/^)  -  <  f(t),^’^>Q  +  g( t)i>>’^(0),  t  >  0,  c 

N  N 

(4.2)  v”(0)  =  P”4. 

where  v^(t)  €  V^,  t  >  0. 

The  state  approximation  given  by  (4,1),  (4.2)  is  the  first  stage  of  our 
approach.  We  argue  that  sufficiently  smooth  solutions  to  (3.6),  (3.7)  are 
approximated  by  solutions  to  (4,1),  (4.2)  with  a  certain  degree  of  uniformity 
in  q.  The  standard  arguments  (see  [31)  yield  the  convergence  of  v^(t)  to  v(t) 
in  the  norm  for  each  t  c  [0,T].  However,  since  the  least  squares 
performance  index  (3,5)  of  particular  interest  to  us  here  involves 
observations  which  are  pointwise  in  the  spatial  variable,  we  will  require  that 
state  approximations  converge  in  the  stronger  norm. 

Theorea  4.1  Suppose  {q  }„  ,  is  a  sequence  in  0  with  q  ♦q^QasN-*-*. 

■■  '  -  — 

Let  v^  denote  the  solution  to  (4.1),  (4.2)  corresponding  to  q^  and  v  the 
solution  to  (3.6),  (3.7)  corresponding  to  q.  Then  if  Hypotheses  (HI)  -  (H5) 
hold,  €  H^(0,1)  and  v  i  H^(0,T;H^(0,1))  we  have 

11m  |v^(t)  -  v(t)|.  *  0 
N-H» 


for  each  t  €  [0,TJ 
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Proof. 

N  N  N  N  N 

Let  q  =  (qpq2,A  )  and  q  -  (qj,q2.^»^)* 

Now 

|v”(t)  -  v(t)|^  <  |v^t)  -  pVt)|j  +  |(P”-I)  v(t)|^. 

The  regularity  assumptions  on  v  and  Hypothesis  (H5)  imply  that  the  second  term 
on  the  right  hand  side  of  the  above  estimate  tends  to  zero  as  N  *  for  each 
t  ^  10, T].  We  therefore  need  only  to  consider  the  term  |v^(t)  -  P^v(t)lj. 
Letting  z^(t)  =  v*^(t)  -  p\(t),  (3.6),  (3.7)  and  (4.1),  (4.2)  yield 

(4.3)  <  z*^(t),)|>^  >Q  +  L(t;q*^)(z^(t),\p^)  *  <  (I-P^)v(t),  >q  +  L  (t;q)(v( t) ,t|»^) 

-  i.(t;q'*)(pVt),4»”),  t  >  0,  /  € 

(4.4)  z”(0)  =  0. 


from  the  above  expression  to  obtain 


Id  /,/  N_  N,2  N,  N,2  ^  .N,  N,2i 


-  v+1  |N|2.v,N,2.,v.  V  .,\ir,N|2 

^T"  *1  2  >0  '0 


J.  1  I/T  I.N^•,2  J.  I  |SN|2  ^  1  i;N,2 

+  ^  Id  -  P  )v|(,  l»ilo  +  2  lij'o 


^  1  iaNi2  .  1  f;N,2  .  d  ,  N  .  N  ^  N  . 
4cj  '^x'o  2  ^*2^  dt  ’^2  ^  ^0 


+  <ApD2Sq  +  <^2*®^%  ■  1^2^^*”^!^* 


.N  ^  N, 


N.  .  N. 


Integrating  (4.6)  from  0  to  t,  applying  Hypothesia  (HI)  -  (H3)  and  recalling 
(4.4)  we  find 


t  T 

0^(t)  <  K  /  9^(s)ds  +  /  p^(s)ds  +  o^(t)  -  o^(0) 
0  0 


where 


9^t)  *  lz”(t)j5  +  |z^t)|^  +  |Dz”(t)lJ, 


N/.x  1  ,/T  oNx  ‘/.mZ  ^  1  ,sN,^.,2  ^  1,  SN,^.,2 
P  (t)  -  ^  Id  -  P  )  v(t)|Q  +  j  l&i(t)lQ  +  jl 


0^(t)  *  q2(t)  <  z^(t),Dz^(t)  >Q  +  <A^(t),Dz’^(t)  >q 


+  <A!?(t),Dz”(t)  >„  -  lA5(t)J,  lz’^(t)J, 


2 
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Once  again  employing  (4.5)  we  obtain 


-N  ^  f  1  ,  N,2  i,x_N,2i  ^  1  ,,N,2 

I"  'o  ^  '=2'’^  'of  '^'o 


j.  ir.  N,2  ^  1  ,.N,2  ^  N,2 

+  C2  |Dz  1q  +  IA^Iq  +  C2  |Dz  Iq 

X  1  fAN,2  ^  ,  N,2 

4c2  ^2*1  ^^2^*  h 


^  V  ,  N,2  ^  oN  ^  .N 

^  4^  'o'"  "3® 


where  “  {v+2)c2  and 


At)  =  4^  l|Aj(t)|J  +  t  4^t)]]\. 


Choosing  C2  such  that  <  1,  we  obtain 


(4.7)  0'^(t)  <  B'^(t)  +-p^/  9”(s)ds 

0 


where 


6”(t)  =  {/  p’^(s)ds  +  a”(0)  +  A^(t)  +  ^  |z^(t)l^}. 


Applying  the  Gronwall  inequality  to  (4.7),  the  desired  result  will  follow 

N 

once  we  have  shown  that  B  (t)  +  0  as  N  +•  for  each  t  t  10,TJ.  The  assumption 
that  V  c  H^(0,T;H^(0,1))  together  with  standard  estimates  yield 


(4.8)  /  p^(s)d8  <  K^f Iq”  -  +  jq^  -  q2li  +  I 

®  N  2 

+  K2l(I  -  P  )vlj  , 


Xl^}  Ivlj 


(4.9)  <  K3{|q5'  -  qjlf  +  Iq”  -  qjl.l  |v(t)|J  +  K^|(I  -  v(t)lj 
and 

(4.10)  &^(0)  <  K3{|q^  -  qj^  +  |q^  -  q^l^l  +  K^I(I  - 

where  1  *  1,2, 3, 4  Is  a  constant  which  does  not  depend  on  N  or  q  c  0. 

It  is  immediately  clear  from  (4.8),  (4.9)  and  (4.10)  that  Hypothesis  (HS) 

together  with  the  assumptions  that  v  €  H^(0,T;  H^(0,1))  and  <>  €  H^(0,1) 

T 

imply  /  p  (s)ds  -►  0,  A  (0)  -►  0  and  A  (t)  -►  0  for  each  t  €  (0,TJ  as  N  +  «. 

0 

Arguments  similar  in  spirit  to  those  above,  although  greatly  simplified,  can 

N  2 

be  used  to  show  (see  [3])  |z  (t)|Q  +0  as  N  +  •*  for  each  t  t  10,TJ.  We 

N 

conclude  therefore  that  6(t)+0asN+<»  for  each  t  €  10, TJ  and  the  theorem 
is  proven. 

An  application  of  the  Sobolev  embedding  theorem  (see  121])  yields  the 
following  corollary. 

Corollary  4.1.  Under  the  hypotheses  of  Theorem  4. 1  we  have 

11m  {  SUP  |v^(t)  -  v(t)l}  ■  0 
N-H»  (0,lj 

for  each  t  [0,T]. 

We  consider  next  the  second  stage  of  approximation;  the  discretization  of 
the  set  of  admissible  parameters,  Q. 

For  each  m~  1,2,...  we  consider  approximation  spaces 
S  C  Cl0,li,  S  *  span  and  S  c  h  (0,T),  S  ■  span  (y  and  let 
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1  12 

ImfJm  denote  mappings  I  :  ClO,ll  +  S  ,  J  :  H  (0,T)  S  .  For  each 
■“  m  m  mm  m 

+  12 

M  *  (m,  ,m„)  with  m.  t  Z  ,  1  *  1,2  define  Q„<^-S  xS  xAxLcrt  hv 

1  i  1  M  —  mj  m2 

(4.11)  Q  =  (Q)  H  I  (Q  )  X  J  (0,)  X  A  X  L. 

n  M  m^  1  ra2  ^ 

We  shall  require  the  following  hypotheses  on  the  mappings  I  : 

M 

(H6)  The  mapping  5  Q  Q.  is  continuous 

ri 

(H7)  For  each  q  t  Q,  (q)  -►  q  as  mj,m2  -*•  •  with  the  convergence 
uniform  in  q  c  Q. 

Note  that  Q  compact  together  with  Hypothesis  (H6)  Imply  that  Ow  Is  a  comnact 

n 

subset  of  0  .  Note  also  that  we  do  not  require  Q,,  c  Q, 

— —  ^ 

1  2 

Once  again  typical  choices  for  S  and  S  are  spaces  of  linear  or  cubic 

m  m 

spline  functions  corresponding  to  meshes  A  *  1/m  and  A  »  T/m  respectively 
with  the  mappings  and  being  the  usual  interpolation  operators*  Under 
appropriate  assumptions  on  Q,  it  is  not  difficult  to  show  that  these  choices 
lead  to  discretizations  which  satisfy  hypotheses  (H6)  and  (H7).  (See  (4), 
151). 

We  consider  the  system  (4.1),  (4.2)  for  q„  ^  0„  .  We  obtain  the  initial 
value  problem  in  R  given  by 

(4.12)  M”w”(t)+  ij  (t;o)w’*(t)  -  F’^(t) 


(4.13) 


w  (0)  -  (  M  )  Wq 


where  a»(a,a,a,a)€fi„ 

M 

and 


a  compact  subset  ofR  xR  xRxR 


(4.14) 


(4.15) 


1  1  x.k  ...N  _.N. 


-  E 
k 


"2  2  k 


.N  ,N. 


N.. v.N, 


3  ..N  .N, 


/  \  x;;;^(t){<D*".^';>^  -  ♦jd )♦■;(!)>  ^  a- 


(4.16)  F”(t)  -  <f(t),*j>Q  +  R(t)*J(l) 


(4.17)  -  <♦,♦';>„ 


l,j  *  l,2,...k^  and 


N 


N,  r  N.  ..N 

V  (t)  »  E  w.(t)^. 

i-1 


We  define  the  sequence  of  approximating  identification  problems  as; 


N  *  N 

(  Pj^)  Find  which  minimizes  over  0^  the  functional  J(q;  v  )  Riven 


by  (3.5)  where  v^  is  the  solution  to  (4.1),  (4.2)  correspondinR  to  q 


or  equivalently: 


Find  a  €  which  minimizes 


N 


J  (o;w’^)  -  E  |Z(5J  -  i  w’J(T(C.;a^  a  ))  ♦?(l)l 


1 
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where  is  the  solution  to  (4.12),  (4.13)  correspondinj?  to  a. 


It  is  Immediately  clear  that  for  each  N,  M  and  a  ^  the  sytem  (4.12), 

(4.13)  admits  a  unique  solution  which  depends  continuously  upon  a.  We  have 

N 

therefore,  that  for  each  N  and  M,  problem  (  )  has  a  solution 

M 

M  *  _M 

(q^)  €  .  Since  there  exists  ^  0  such  that 


N 

(«!«)  “  (q»j)  •  Now  Q  is  compact.  Therefore,  there  exist  {q„  ) 
n  M  n  M, 

N  ^ 

*  — r  4  it 

and  an  element  q  t  Q  such  that  q,,"^  q  in  (2  as  i  ,k  +  «.  This 

\ 

subsequence  can  always  be  chosen  so  that  +  •  and  ' 

j  ,k  -►  *.  Now 


/~N  •! 


N.  ^  N 


1..  ^  i,,  1., 

J((qf,j^)  ;v  ')  <  J(q;v  '),  q  ^  0^^ 


and  consequently,  from  (4.11)  ,  we  have 


(4.18)  J((qj^)  ;v  ^)  <  J(  (q);v  '),  q  c  Q. 


Hypothesis  (H7)  and 


N.  N, 


l<X>  -  "  is  ‘I  *  lx  -  ’  'o 


i  *  * 

imply  (qjX)  q  in  (i  as  .1,k  ♦  •.  Taklnq  the  limit  as  i,k  •  in  (4.18) 

K 

and  applying  Corollary  4.1  (with  an  appropriate  re-lndexini?  and  the  assumption 
that  the  necessary  rei^ularity  conditions  are  satisfied)  we  obtain 
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*  * 

J(q  ;v  )  <  J(q;v),  q  €  Q 

where  v*  denotes  the  solution  to  (3.6),  (3.7)  corresponding  to  q*;  or  that  q* 
is  a  solution  to  problem  (  P  ). 

We  summarize  these  results  in  the  following  theorem. 


N 


Theorem  4.2  For  each  N  »  1.2....  and  M  €  x  let  problem  (  P.,  )  be  as 


it  has  been  defined  above.  Then  if  hypotheses  (HI)  -  (H?)  hold,  each  problem 
N 

N  *  N  * 

(  Pj^)  has  a  solution  (dj^)  •  The  sequence  f(d^)  )  admits  a  ^-convergent 


subsequence,  the  limit  of  which,  q*,  is  an  element  of  Q.  If  c  h\o,1)  and 


V*,  the  solution  to  (3.6),  (3.7)  corresponding  to  q*  is  an  element  in  H^(0,T; 


H^(0,1))  then  q  is  a  solution  to  problem  (P  ).  Moreover,  for  anv  convergent 

*  N  *  *  -* 

subsequence  Hq^/)  )  ^  ^  ♦  d  .  +  •  and  (rap^^, 


— 4r  1  I  ifr 

(n>2^k  *  as  j,k  •  and  v  €  H  (0,T;H  (0,1))  we  have  that  q  is  a  solution 


to  problem  ( P  ) 


5.  Nuaerical  Results 


In  this  section  we  present  a  brief  summary  of  some  preliminary  numerical 
results.  Our  primary  Intent  here  Is  to  simply  demonstrate  the  feasibility  and 
efficacy  of  our  approach.  A  more  complete  study  involving  a  somewhat  broader 
spectrum  of  examples  (Including  pathologies)  and  more  sophisticated  approaches 


to  solving  the  approximating  optimization  problems  will  be  discussed  in  detail 
elsewhere. 

We  consider  inverse  problems  for  systems  of  the  form  (2.3)-{2.6) 
involving  the  estimation  of  a  spatially  varying  diffuslvity  coefficient 
P=  P(x).  We  assume  that  the  advection  rate  V  Is  constant  in  time,  that  there 
is  no  boundary  flux  and  that  all  parameters,  with  the  exception  of  the 
diffuslvity,  are  known.  We  also  allow  for  the  inclusion  of  a  sink/source 
term,  F  =  F(t,x)  in  (2.3). 


Transforming  both  the  space  and  time  coordinates  to  dimensionless 
variables;  y  =  x/A,  s  -  Vt/i,  we  obtain 


(5.1)  v(s,y)  =  {q(y)  v(s,y)}  ~  v(s,y)  -  a  v(s,v)  +  f(s,v),  s  >  0,  vt  (0,1' 


(5.2)  -q(0)  I-  v(s,0)  +  v(s,0)  »  0 


(5.3)  -q(l)  ~  v(s,l)  -  0 


(5.4)  v(0,y)  -  4»(y) 


V  ^  (0,1) 


where  v(s,y)  *  u(ls/l',fy),  q(y)  -  P<ty)/<{/,  a  *  f(s,v)  -  «F(£s/ l',«v)/l/. 
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and  <Ky)  =  ♦(J.y). 

In  addition,  for  our  discussions  here,  we  assume  that  observations 

{C..}.  ,  for  V  (equivalently,  u)  at  times  s^,  1  =  and 

. 

positions  y.  j  >  l,...,v  have  been  provided  and  consider  the  least  squares 

J  9 

performance  Index  of  the  form 

W  V 

J(q;v)  -  I  I  |v(s.,y,)  -  e..i 
1-1  1-1  ^ 

where  v  is  the  solution  to  (5.1)  -  (5.4)  corresponding  to  q. 

All  computations  were  performed  on  IBM  3081  processors  at  Brown 

University  and  the  University  of  Southern  California.  The  approximating  state 
N 

spaces  V  were  chosen  as  the  span  of  cubic  B-splines  defined  with  respect  to 

the  uniform  partition  {0,1/N,  2/N,...,l}  of  the  interval  t0,ll.  The  set  of 

admissible  parameters  was  discretized  using  linear  spline  functions  with 

respect  to  the  mesh  {0, 1 /M,2/M, . . . , 1 }.  Note  that  in  this  case  we  have 

=  N  +  3  and  =  M  +  1.  The  Inner  products  In  (4.14)  -  (4.17)  were 

computed  using  a  composite  two  point  Gauss-Legendre  quadrature  formula.  The 

use  of  B-spllne  bases  leads  to  banded  generalized  mass  (M  )  and 
N 

stiffness  (/.„)  matrices. 

n 

The  finite  dimensional  optimization  problems  (P^)  were  solved  using  an 
Iterative  Levenberg-Marquardt  scheme  as  implemented  in  the  IMSL  Library 
routing  ZXSSQ.  Gradients  and  Jacobians  are  computed  bv  the  routine  using 
finite  difference  approximation  with  rank  one  updates  in  each  iteration.  We 
also  attempted  to  solve  the  finite  dimensional  optimization  problems  using  a 
quasi-Newton  algorithm,  a  scheme  due  to  Broyden,  Fletcher,  Goldfarb  and  Shanno 
(BFGS)  (see  19]),  with  analytical  gradients  computed  using  a  co-state 
formulation  (see  [2],  [8]).  Our  preliminary  findings  point  to  the  conclusion 
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that  the  latter  approach  Is  Inferior.  We  found  that  It  was  extremely 
difficult  to  obtain  accurate  search  directions  with  the  gradients  computed  in 
this  manner.  However,  we  are  continuing  to  investigate  these  ideas  with  the 
Intent  of  improving  numerical  performance. 

In  the  examples  that  follow,  the  compactness  constraints  on  the  set  of 
admissible  parameters  Q  were  not  explicitly  enforced  when  the  finite 
dimensional  optimization  problems  were  solved.  As  could  be  expected,  this  did 
lead  to  some  conditioning  problems  when  M  became  large.  The  use  of  a 
constrained  optimization  package  to  solve  the  finite  dimensional  optimization 
problems  is  currently  under  study. 

k"^ 

The  Initial  value  problems  (4.12),  (4.13)  in  R  were  solved  in  each 
iteration  using  Gear's  method  (IMSL  Library  routine  OGEAR)  for  stiff  systems. 

EXAMPLE  5.1 

We  consider  the  system  (5.1)  -  (5.4)  with 

q(y)  *  1  +  y^, 

f(s,y)  -  (2  -  8y  +  /y^)e  ®, 

♦(y)  -  2  +  2y  -  y^ 

and  0*0.  We  estimated  q  from  observations  at  {(s  ,  V,)},  «  i  o  o 
generated  using  the  true  solution  to  the  system  2,...  ,4 

(5.5)  v(8,y)  *  (2  +  2y  -  y^)e  ® 

where  Sj^  *  .251  and  y^  -  .25j.  The  initial  estimate  for  q  supplied  to  the 


optimization  routine  was  taken  to  be 


q  (y)  =1,  0  <  y  <  1. 


Our  results  with  N  =  32  and  various  values  of  M  are  summarized  in  Table  5.1 
below.  The  case  M  =  0  corresponds  to  the  best  fit  taking  q  to  be  constant 
over  the  entire  interval. 


* 

* 

* 

* 

* 

qo(y) 

qj(y) 

q3(y) 

q^(v) 

q(y) 

1.172 

0.961 

0.978 

0.987 

1.172 

0.992 

1.006 

1.014 

1.015 

1.010 

1.172 

1.065 

1.052 

1.050 

1.042 

1.040 

1.172 

1.097 

1.085 

1.092 

1.090 

1.172 

1.209 

1.143 

1.164 

1. 164 

1.160 

1.172 

1.208 

1.264 

1.236 

1.250 

1.172 

1.354 

1.354 

1.364 

1.363 

1.360 

1.172 

1.426 

1.500 

1.479 

1.490 

1.490 

1.172 

1.498 

1.645 

1.624 

1.637 

1.640 

1.172 

1.571 

1.791 

1.769 

1.804 

1.810 

1.172 

1.643 

1.937 

1.914 

1.972 

-3 

„  -4 

-6 

-6 

6.  X  10 

2.  X  10 

4.  X  10 

3.  X  10 

• 

X 

o 

< 

0:49.92 

1:00.53 

1:29.49 

2:02.02 

2:35.82 

TABLB  5.1 


EXAMPLE  5.2 


In  this  example  we  estimate  the  non-monotone,  single  neak  function  q 
given  by 


q(y)  *  (1  +  ay)e 


-by 
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with  a  “  10.43  and  b  *  2.43.  Takln*;  a  ■  0, 


f(s,y)  ■  -  2ab)  +  y(2ab  +  4a  —  2b)  +  (2b  -  2a  +  2  -  4e  ^))e  \ 


4>(y)  -  2  +  2y  -  y  . 


The  true  solution  to  the  system  Is  again  given  by  (5.5).  Using  the  same 
observation  points  as  In  the  previous  example,  setting 


q^(y)  -  1.5,  0  <  V  <  1 


and  taking  N  =  6  and  M  -  4,  we  obtained  the  fit  for  q  shown  In  Figure  5.1 
below. 


2.00 


(q!)* 


1.00* - ' - > - 1 - - 

0.0  .25  50  .75  ftOO 


TXGUiB  3*1 


The  residual  was  J((q^)  )  ■  3.  *  lO”^. 


-38- 


The  scheme  performed  well  when  N  and  M  were  chosen  relatively  small. 
Increasing  either  N  or  M  Independently  yielded  some  Imorovement  Initially. 
However,  this  was  followed  by  the  onset  of  ill  conditioning  in  the  finite 
dimensional  optimization  problems  as  the  number  of  degrees  of  freedom  in 
either  the  state  or  the  parameter  space  discretizations  were  increased. 

(Based  on  our  initial  computational  findings  with  other  examples,  we 
anticipate  that  enforcing  the  compactness  constraints  on  the  set  of  admissible 
parameters  will  remedy  this  situation.)  The  scheme  also  began  to  have  some 
difficulty  if  the  parameters  a  and  b  were  chosen  so  as  to  cause  the  total 
variation  of  q  to  be  too  large. 

EXAMPLE  5.3 

We  estimate  a  depth  dependent  bioturbation  coefficient  using  a  set  of 
observations  from  a  volcanic  ash  concentration  profile  measured  in  a  sediment 
core  sample  taken  in  the  North  Atlantic.  We  consider  the  model  given  by  (5.1) 
-  (5. A)  with  a  ■  0  (we  are  dealing  with  a  non-radioactive  or  conservative 
tracer),  f  ■  0  and 

(5.6)  ♦(y)  -  md(y) 

where  6  denotes  the  unit  impulse  at  zero  and  m  is  the  total  mass  of  ash  in  the 
sample  as  determined  from  the  data. 

We  were  provided  with  the  estimate  V-  2.5  cra/kyr  for  the  sedimentation 
rate  in  Che  region  where  the  core  sample  from  which  our  ash  data  came  was 
taken.  In  addition,  by  invoking  (possibly  inappropriately  in  the  presence  of 
depth  dependent  mixing)  the  observation  in  [llj  Chat  the  concentration- 
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welghted  mean  depth  of  the  ash  profile  In  the  historical  layers  is  the  depth 
in  the  core  at  which  the  ash  layer  would  have  been  observed  had  no  mixing 
taken  place,  we  estimate  1,  the  depth  to  which  mixing  takes  place  (see  {6])  to 
be  17.23  cm. 

Using  the  estimates  for  V  and  I  given  above  we  are  able  to  convert  the 
depth  scale  on  which  our  data  is  specified  to  an  equivalent  s-scale.  Our 
observations  turn  out  to  be  given  at  temporal  locations  Sj^  *  .1449251, 

1  *  0,1, 2,. ..,16  and  of  course  all  at  the  spatial  location  y  ■  1. 

We  set  N  >  32  and  used  our  scheme  to  estimate  q  for  various  values  of 
M.  We  approximated  the  Impulse  initial  conditions  (5.6)  by 

4»(y)  -  m^2^(y) 

^32  ^  ^32 

where  41,  denotes  the  normalized  (/  ♦<,  *  1)  cubic  B-spllne  corresponding  to 

the  uniform  mesh  {0,1/32,2/32,...!}  which  is  centered  over  zero.  This 
approximation  is  {ustlfled  by  the  relatively  narrow  support  of  the  B-snlines 
and  the  fact  that  it  eliminates  the  error  which  would  be  Introduced  if  anv 
other  impulse  function  approximation  were  proiected  onto  the  subsoace  of 
splines. 

The  initial  estimate  for  q  was  taken  to  be 

q°(y)  -  .02198,  0  <  y  <  1, 

the  optimal  estimate  for  q  obtained  in  [6]  using  this  set  of  observations  and 
the  assumption  that  the  mixing  intensity  is  constant  throughout  the  mixed 
layer. 


The  model  (3.1)  -  (5.4)  with  deoth  dependent  mixing  rate  yielded  some, 
although  not  a  significant,  reduction  In  the  residual  when  compared  to  the  fit 
obtained  using  a  constant  mixing  rate  model.  (See  Figure  5.2  below.) 
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In  addition  the  optimal  mixing  rate  profiles  produced  by  our  scheme  did  not 
agree  with  the  widely  accepted  hypothesis  that  mixing  Is  most  Intense  near  the 
sea  floor/sea  Interface  and  then  decreases  with  depth  (see  Figure  5.3). 


^  V 
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FIGDIE  5.3 

Based  upon  these  findings,  we  conclude  therefore  that  the  inclusion  of  a 
depth  dependent  mixing  rate  alone  can  not  significantly  improve  our  ability  to 
explain  this  set  of  observations.  Other  enhancements  of  the  original  model, 
(e.g.,  time  dependent  sedimentation  rate,  porosity  and/or  compact! fi cat ion 
efffects)  must  be  considered.  The  scheme  developed  here  should  prove  to  be  a 
valuable  tool  in  the  investigation  and  evaluation  of  these  modeling  ideas. 


ACKMUHLEDGEMENT  The  authors  would  like  to  gratefully  acknowledge 
Nr.  Chunming  Wang  for  the  assistance  he  provided  in  carrying  out  the 
computations  discussed  in  this  paper  and  Professors  John  Imbrie  and  Warren 
Prell  of  the  Department  of  Geological  Sciences  of  Brown  University  for  their 
continued  encouragement  and  stimulating  remarks  on  this  work.  They  would  also 
like  to  thank  them  for  supplying  the  volcanic  ash  data  from  the  laboratory  of 
Dr.  William  Ruddiman  of  The  Lamont-Doherty  Geological  Observatory  of  Columbia 
University  which  was  used  in  Example  5.3. 


s.  .. 


■EraRENCES 


[IJ  H.T.  Banks,  On  a  variational  approach  to  some  parameter  estimation 

problems,  ICASR  Report  No.  85-32,  Institute  for  Computer  Applications  in 
Science  and  Engineering,  NASA  Langley  Research  Center  Hampton,  VA,  June 
1985. 

12J  H.T.  Banks  and  J.M.  Crowley,  Estimation  of  material  parameters  In 

elastic  systems,  LCDS  Report  No.  84-20,  Brown  University,  June  1984. 

[3]  H.T.  Banks,  P.  Karelva  and  P.D.  Lamm,  Estimation  techniques  for 
transport  equations,  in  Mathematics  in  Biology  and  Medicine 
(Proceedings,  Bari  1983),  V.  Capasso,  et.  al. ,  eds.,  LN  in  Blomath  Vol. 
57,  Springer,  New  York,  1985,  428-438. 

[4]  H.T.  Banks,  and  P.D.  Lamm,  Estimation  of  variable  coefficients  in 
parabolic  distributed  systems,  IEEE  Trans.  Auto.  Control,  AC-30,  (1985), 
386-398. 

15]  H.T.  Banks  and  K.A.  Murphy,  Estimation  of  coefficients  and  boundary 

parameters  In  hyperbolic  systems,  LCDS  Technical  Rept.  No.  84-5,  Brown 
University,  February,  1984;  SIAM  J.  Control  and  Opt.,  to  appear. 

[6J  H.T.  Banks  and  I.G.  Rosen,  Fully  discrete  approximation  methods  for  the 
estimation  of  parabolic  systems  and  boundary  parameters,  LCDS  Tech.  Rep. 
No.  84-19,  Brown  University,  May,  1984,  Acta.  Applic.  Math.,  to  appear. 

17J  L.K.  Benninger,  R.C.  Aller,  J.K.  Cochran  and  K.K.  Turekian,  Effects  of 
biological  sediment  mixing  on  the  210Pb  chronology  and  trace  metal 
distribution  In  a  Long  Island  Sound  sediment  core.  Earth  Planet.  Scl. 
Lett.  ^(1979),  241-259. 

(8J  G.  Chavent,  Identification  of  distributed  parameter  systems:  About  the 
output  least  squares  method.  Its  implementation,  and  Identif lability, 
Proc.  5th  IFAC  Symp.  on  Identification  and  Systems  Parameter  Estimation. 
Pergamon  Press,  1979,  85-97. 

I9J  R.  Fletcher,  Practical  Methods  of  Optimization  John  Wiley,  New  York, 
1981. 

IIOJ  A.  Friedman,  Partial  Differential  Equations  of  Parabolic  Type.  Prentice- 
Hall,  Englewood  Cliffs,  1964. 

[11]  N.L.  Gulnasso,  Jr.,  and  D.R.  Schink,  Quantitative  estimates  of 
biological  mixing  rates  in  abyssal  sediments,  J*  Geophys.  Res.  80(1975), 
3032-3043. 

[12]  D.  Kadko  and  G.R.  Heath,  Models  of  depth-dependent  bloturbatlon  at  MANOP 
site  H  in  the  eastern  equatorial  Pacific,  J.  Geophys.  Res.  89(1984), 
6567-6574. 


-43- 


[13]  R.S.  Keir,  Recent  Increase  in  Pacific  CaCO^  dissoultlon:  a  mechanism 
for  generating;  old  ages.  Mar,  Geol,  59(19B4) ,  227-250. 

1 14]  J.L.  Lions,  Optimal  Control  of  Systems  Governed  by  Partial  Differential 
Equations,  Springer  Verlag,  New  York,  1971, 

[15]  D.R.  Lynch  and  C.B.  Officer,  Nonlinear  Parameter  estimation  for  sediment 
cores,  Chera.  Geol.  44(1984 ) ,  203-225. 

[16]  C.B,  Officer  and  D.R.  Lynch,  Interpolation  procedure  for  the 
determination  of  sediment  parameter  from  time-dependent  flux  incuts, 
Earth  Planet.  Scl.  Lett.  61(1982).  55-62. 

[17]  T.-H.  Peng  and  W.S.  Broecker,  The  impacts  of  bloturbatlon  on  the  aee 
difference  between  benthic  and  planktonic  foraminifera  in  deep  sea 
sediments,  Nucl.  Instr,  and  Meth.  B5  North  Holland,  1984,  346-352. 

[18]  T,-H.  Peng,  W.S.  Broecker  and  W.H,  Berger,  Rates  of  benthic  mixing  in 
deep-sea  sediment  as  determined  by  radioactive  tracers,  Ouaternary  Res. 
n_(1979),  141-149. 

[19]  W.F.  Ruddlman  and  L.K.  Glover,  Nixing  of  volcanic  ash  zones  in  subpolar 
North  Atlantic  sediments,  in  The  Ocean  Floor  (R.A.  Scrutton  and  M. 
Talwanl,  eds),  John  Wiley,  1982,  37-60, 

[20]  M.H.  Schultz,  Spline  Analysis,  Prentice  Hall,  Englewood  Cliffs,  1973. 

[21 J  R.G.  Showalter,  Hilbert  Space  Methods  for  Partial  Differential 
Equations.  Pitman,  London,  1977, 

[22]  B.K.  Swartz  and  R.S,  Varga,  Error  bounds  for  Spline  and  L-Spline 
interpolation,  J.  Approx.  Theory,  M1972),  6-49, 

[23]  H.  Tanabe,  Equations  of  Evolution.  Pitman,  London,  1979. 

[24]  K.K.  Tureklan,  J.K.  Cochran,  and  D.J.  DeMaster,  Bloturbatlon  in  deep-sea 
deposits:  Rates  and  consequences,  Oceanus  21(1978),  34-41. 


